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Abstract 

In a previous work we developed a convex infinite dimensional linear program- 
ming (LP) approach to approximating the region of attraction (ROA) of polynomial 
dynamical systems subject to compact basic semialgebraic state constraints. Finite 
dimensional relaxations to the infinite-dimensional LP lead to a truncated moment 
problem in the primal and a polynomial sum-of-squares problem in the dual. This 
primal-dual linear matrix inequality (LMI) problem can be solved numerically with 
standard semidefmite programming solvers, producing a hierarchy of outer (i.e. ex- 
terior) approximations of the ROA by polynomial sublevel sets, with a guarantee 
of almost uniform and set-wise convergence. In this companion paper, we show 
that our approach is flexible enough to be modified so as to generate a hierarchy of 
polynomial inner (i.e. interior) approximations of the ROA with similar convergence 
guarantees. 

1 Introduction 

Given an autonomous nonlinear system and a target set, the region of attraction (ROA) 
is the set of all states that end in the target set at a given time without leaving the state 
constraint setj^j The ROA is one of the principal sets associated to any dynamical system 
and goes by many other names in the literature (e.g., backward reachable set or capture 
basin [4J). 

In [6j we showed (in a controlled setting) that there is a genuinely primal convex character- 
ization of the ROA. Optimization over system trajectories is formulated as optimization 
over occupation measures, leading to an infinite dimensional linear programming (LP) 
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problem in the cone of nonnegative measures. Finite dimensional relaxations of the dual 
of this problem then provide a converging sequence of outer approximations to the ROA. 
For a description of alternative techniques for numerical approximations of the ROA, 
please consult [5] or [6] and the many references therein. 

In this paper we show, within the same measure-theoretic framework, that there ex- 
ists an infinite dimensional LP whose finite-dimensional relaxations provide a converging 
sequence of inner approximations to the ROA. This paper can therefore be seen as a com- 
plement to [6j. To simplify our developments and to emphasize our contribution, we focus 
only on the uncontrolled setting. The main idea is to construct a converging sequence of 
outer approximations to the complement of the ROA. There are certain difficulties, topo- 
logical in nature, associated with this approach. A careful distinction had to be made 
between trajectories leaving the constraint set and trajectories hitting its boundary. This 
then translates to a (sometimes subtle, but necessary) distinction between open and closed 
semialgebraic sets. Fortunately, the LP formulation proposed in [6j was flexible enough 
to allow for these modifications. 

Generally speaking, and consistently with our previous work [6j, we believe that the main 
virtues of our approach are overall convexity, conceptual simplicity and compactness. Both 
primal and dual finite-dimensional relaxations turn out to be linear matrix inequalities 
(LMI), also called semidefinite programming (SDP) problems, with no tuning parameters 
besides the relaxation order and no initialization data besides the defining ingredients of 
the problem. In addition, the inner approximations obtained are particularly simple 
they are given by a sublevel set of a single polynomial of a predefined degree. Therefore, 
an ROA approximation in analytic form can be readily obtained by solving a single LMI 
using freely available software (e.g., SeDuMi [T2]). 



2 Problem statement 

Consider the autonomous system 

x(t) = f(t,x(t)), x{t)eXcR n , te[0,T] (l) 

with a given vector field / with polynomial entries fa G R[t, x], i = 1, . . . , n, final time 
T > 0. The state trajectory x(-) is constrained to a nonempty operj^] basic semialgebraic 
setH 

X :={xeR n : g x (x) > 0}, (2) 
where the polynomial g x £ R[x] is such that the set 

X :={xeR n : g x {x) > 0} D X 

is compact^) 

6 The requirement of the constraint set being open is merely technical, for this considerably simplifies 
the developments and the proofs. 

7 For clarity of exposition we consider the constraint set given by a single superlevel set of a polynomial. 
The approach can, however, be straightforwardly extended to constraint sets defined by the intersection 
of finitely many polynomial superlevel sets. 

8 Note that the closed semialgebraic set X — {x : gx{x) > 0} can be strictly larger than the closure 
of the open semialgebraic set X — {x : Qx{ x ) > 0} 5 consider in K. e.g. Qx{ x ) — (1 — x 2 )(2 + x) 2 . For 
a similar reason, note also that X bounded does not imply X bounded. Indeed, in M 2 with gx(x) — 
(1 — x\ — #2) (2 + xi) 2 we have X — {x : ||x|| < 1} and X — {x : ||x|| < 1} U {x : X\ = -2}. 
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The vector field / is polynomial and therefore Lipschitz on the compact set X. As a 
result, for any xq G X there exists a unique maximal solution x(-) to ODE ([!]). The time 
interval on which this solution is defined contains the time interval on which x(t) G X. 

2.1 Region of attraction (ROA) 

Given a final time T and an open bounded basic semialgebraic target set 

X T := {xeW 1 : g T (x) > 0} C X, 
the region of attraction (ROA) is defined as 

Xq := jx G X : 3x{-) s.t. x(t) = f(t,x(t)), x(0) = x Q , x{T) G X T ,x(t) G X,V£ G [0,T]J 

(3) 

In words, the ROA is the set of all initial states from X for which the unique solution 
to stays in X for all t G [0, X] and ends in the target set at time T. 

2.2 Complement ROA 

The idea to get inner approximations of the ROA X is to construct outer approximations 
of the complement ROA XJ := X \ X . By continuity of solutions to Q, the set Xq is 
equal to 

= [x G X : 3x(-) s.t. x(t) = f(t,x(t)) and 

3t G [0,T] s.t. x(t) G X a and/or x(T) G X£}, 

where 

XJ := {xeR n : > 0, g T (x) < 0} 

is the complement of Xt in X and 

X a := {x G R n : g x {x) = 0}. 

In words, Xq is the set of initial states that give rise to trajectories which do not end up 
in Xt at time T and/or violate the state constraint at some point between and T. 

3 Occupation measures 

In this section we introduce the concept of occupation measures and show how the non- 
linear system dynamics can be equivalent ly described by a linear equation on measures. 

Notation We will use the following notation. The vector space of all signed Borel 
measures with support contained in a Borel set K is denoted by M(K). The support (i.e., 
the smallest closed set whose complement has a zero measure) of a measure fi is denoted 
by spt fi. The space of continuous functions on K is denoted by C{K) and likewise the 
space of continuously differentiable functions is C l {K). The indicator function of a set K 
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(i.e., the function equal to one on K and zero otherwise) is denoted by Ik{')- The symbol 
A denotes the n-dimensional Lebesgue measure (i.e., the standard n-dimensional volume). 
The integral of a function v w.r.t a measure ji over a set K is denoted by J R v(x) dfi(x). 
Sometimes we for simplicity omit the integration variable and/or the set over which we 
integrate if they are obvious from the context. 

Now assume x G X and define the first hitting time of Xq as 

r(x ) := min{T,inf{£ > : x(t \ x ) G X d }}, (4) 

where x(- \ x ) denotes the unique trajectory starting from x (which is well defined on 
the time interval [0,r(xo)]). Then we define the occupation measure associated to the 
trajectory starting from xq by 

t(x ) 

dt 



H(A X B | Xq) := / I A xB{t,x{t)) 
Jo 



for all Borejjsets A x B C [0, T] x X. The interpretation is that the occupation measure 
measures the time spent by the trajectory x(- \ x ) in subsets of the state space. 

The occupation measure enjoys the following important property: for any measurable 
function v(t,x) the equality 



I v(t,x(t))dt= / v(t, x) d/l(t, X | Xq) (5) 

J0 J[0,T]xX 



holds. In words, the time integral of a function evaluated along the trajectory x(- \ x ) 
is equal to the integral of the function w.r.t. the occupation measure associated to x . 
Therefore, loosely speaking, all information about the trajectory x(- \x ) is encoded by 
the occupation measure fi(- \ x Q ). 

Now suppose that the initial state is not a single point but that its spatial distribution is 
given by an initial measure fi G M(X). Then we define the average occupation measure 
[i G M([0,T] x X) as 

fi(A x B) := / fi(A x B \ x ) dfi (x ). 

J X 

Lastly, we define the final measure fir £ M([0,T] x X) by 

fi T (B) := / I B {x(T | x ))dfio(x ). 
Jx 

To derive an equation linking the three principal measures, consider a test function v G 
C 1 ([0, T] x X) evaluated along a trajectory. Using the chain rule and equation (0) we 



} For brevity we drop the adjective "Borel" in the sequel. 
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obtain 



rr(x ) ^ 

'(t(x ), x(t(x ) I x )) - v(0, x ) = J ~J t v & x ( f I dt 



'o 



— + grad v ■ f(t, x(t \ x )) ) dt 



[0,T]xX 



dv \ 
— + grad v • f(t,x) J dfi(t, x \ x ) 

Cv(t, x) d/i(t, x | xq) 

[0,T]xX 

where the linear operator C : C l (^),T} x X) —> C([0,T] x X) is defined by 

v i — y JL,v '.= — — h gradv • /. 
at 

Integrating the above equation w.r.t. fi leads to the equation 
/ v(t, x) dfiT{t, x) — \ i;(0, x) dfio(x) — / Cv(t,x) dfi(t,x) V v G C 1 ([0, T] x X), 

J[0,T]xX JX J[0,T]xX 

(6) 

which is a linear equation linking the measures fi 0l fi and fiT- Equation Q is sometimes 
referred to as Liouville's equation. 



4 Primal LP 

In this section we follow the approach developed in [6j and derive an infinite-dimensional 
linear programming (LP) characterization of the complement ROA Xq. Certain sublevel 
sets of feasible solutions to the dual of this LP then yield inner approximations to the 
ROA X . 

The basic idea is to maximize the mass of the initial measure fio under the constraint 
that it is dominated by the Lebesgue measure, i.e., fio < A. System dynamics is captured 
by Liouville's equation ^ and state and terminal constraints are handled by suitable 
constraints on the support of the measures. The key idea is then to split the final mea- 
sure in two measures such that each measure is supported on a suitable compact basic 
semialgebraic set. More explicitly, we let 

Mt := lA + /4 

with \i\ G M([0,T] x X d ) and /i| G M({T} x X%). That is, we require that spt \i\ C 
[0, T] x Xq and spt/iy C {T} x The interpretation is that measure models the 
trajectories that leave X 1 whereas measure \3p models the trajectories that end in X^ 
(i.e., not in Xt)- These support constraints on the final measure(s) along with system 
dynamics enforce that the support of the initial measure fio must be contained in Xq. 
Since there are no other constraints on /i besides fio < A, maximization of its mass 
should yield the restriction of the Lebesgue measure A to Xq. 

The constraint fio < A can be rewritten equivalently as fi + fio — A for some nonnegative 
slack measure fio G M(X). This is equivalent to requiring that J w dfio + J w dfio — J w dX 
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for all test functions w G C(X). In addition, we can drop the time argument from the 
definition of Ht since its time component is supported on a singleton. 

This leads to the following optimization problem: 

p* — sup J 1 dfiQ 

s.t. fvd^ + J v(T, -)dfil - J v(O r )d/i = JCvd/i VuG C^QC^T] x X) 
J w dfiQ + J w dfiQ — J w dX V w G C (X) 

fi>o > 0, fi> > 0, /iy > 0, > 0, /} > 
spt /i C [0, T] x X, spt no C X, spt fio C X 
spt //y C [0, T] x X 5 , spt /iy C Xy, 

(7) 

where the supremum is over the vector of nonnegative measures 

Oo, Ht, l4, M G m W x M([0,T] x X) x M([0,T] x X a ) x M(X£) x M(X). 

Problem ([7]) is an infinite-dimensional LP in the cone of nonnegative measures. Indeed, 
the objective is linear, the first two constraints are linear equality constraints and the 
remaining constraints are conic constraints (the set of nonnegative measures supported 
on a given set is a positive cone in the vector space of all measures supported on the same 
set). 

The discussion leading to problem is formalized in the following result. 

Theorem 1 The optimal value of LP problem fijty is equal to the volume of the comple- 
ment ROA Xq ; that is, p* — A(Xq). Moreover, the supremum is attained by the restriction 
of the Lebesgue measure to the complement ROA Xq. 

Proof: Closely follows arguments in [6j. By definition of the relaxed complement ROA, 
the unique trajectory x(-) associated to any initial condition xo G Xq either hits Xq at 
some t G [0,T] or ends in X£. Therefore for any initial measure Ho < ^ with spt Ho C 
X C X there exist an occupation measure /i, final measures Hti f^r an< ^ a sl ac k measure 
fio such that the constraints of problem Q are satisfied. One such measure Ho is the 
restriction of the Lebesgue measure to Xq, and therefore p* > A(XJ). 

Now we show that p* < A(Xq). Take a vector of measures (ho, V, Vt, V>h M feasible 
m and suppose that A (spt ho \ Xq) > 0. Since any level set of a polynomial has a zero 
Lebesgue measure we have \(Xg) = and 

A(spt ho \ (X C q U X d )) = A(spt ho \ X%) > 0. 

By a superposition principle [TJ Theorem 3.2] using arguments of [61 Appendix A, Lemma 4], 
there exists a family of admissible trajectories of the ODE ([I]) starting from Ho generat- 
ing the occupation measure h an d the final measure ht = + Mt- However, this is 
a contradiction since spt Ho \ {Xq U Xq) C X , which means that all trajectories start- 
ing from spt Ho \ {Xq U Xq) neither hit Xq nor end in X T . Thus, A(spt/i \ XJ) = 
and so A(spt/i ) < ^PQ))- Combining this with the constraint /i < A we get Ho(X) = 
Ho (spt Ho) < A(spt/io) < ^(^o) fe r an y feasible Ho- Therefore p* < A(Xq) and thus in 
fact p* = A(X C ). □ 
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5 Dual LP 



In this section we derive a dual LP on continuous functions, prove the absence of a 
duality gap between the primal and dual LPs and relate feasible solutions to the dual to 
the indicator function of the complement ROA Xq. 

By standard infinite-dimensional LP theory (see, e.g., [2j), the dual to LP Q reads 

d* = inf j x w(x) d\(x) 
s.t. Cv{t 1 x) < 0, 

w(x) > i>(0, x) + 1, 
v(T,x) > 0, 
v(t,x) > 0, 
w(x) > 0, 

where the infimum is over (v,w) G C 1 ([0,T] x X) x C(X). 

The intuition is that given x G Xq the constraint Cv < forces v to decrease along 
trajectories as long as it does not hit Xq or end in X£. Because of the constraint v > 
on [0, T] x Xq U {T} x Xf we must have t>(0, •) > on Xq. Consequently, w(x) > 1 on 
Xq. This instrumental observation is formalized in the following Lemma. 

Lemma 1 If Cv < on [0, T] x X, v > on ([0, T] x X 5 )U({T} x Xj) andw > v(0,-) + l 
on X , then w > 1 on Xq. 

Proof: Take x G Xq and consider the first hitting time of Xq, t := r(x ), defined 
by Ml). By definition of Xq the trajectory starting from xq will either hit or end in 
X£. Therefore x(r) G ([0,T] x X a ) U {{T} x X^) and x(t) G X for t G [0,r]. Therefore 
u(r,x(r)) > 0, Cv(t,x(t)) < 0, Vt G [0,r] and so 

< v(r, a;(r)) = v(0, xq) + / ^(^)) dt < ^(0, x ) < w(^o) — 1- 

Jo 

□ 

The following result is of key importance for subsequent developments. 

Theorem 2 There is no duality gap between primal LP problems |?p on measures and 
dual LP problem (Q) on functions, that is, p* = d*. 

Proof: Here we only outline the basic steps; for a detailed argument in a similar setting 
see Theorem 2). Since the supports of all measures are compact, the initial measure is 
dominated by the Lebesgue measure and the final time is finite, we have fio(X) < A(X) < 
oo, /i T ([0,T] x X) = no(X) < oo and /i([0,T] x X) < 7> T ([0,T] x X) < oo, where 
the last two inequalities follow by plugging in v(t,x) = 1 and v(t,x) = t in Liouvillel's 
equation Q. Therefore p* < oo and the feasible set of problem Q is weakly-* bounded. 
Furthermore, the feasible set of Q is nonempty since (/iq, fi, /i^, /i^, (lq) — (0, 0, 0, 0, A) is 
a trivial feasible point; therefore < p* < oo. The absence of a duality gap then follows 
from [21 Theorem 3.10] using Alaoglu's theorem (see, e.g., [TQl Chapter 5]) and the weak-* 
continuity of the adjoint of the operator C. □ 



\/(t,x,u) G [0,T] x X 
Vi G X 
\/x G Xff 

V(t, x) G [0, T] x X a 
Vx G X, 



(8) 
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Next, we establish our first convergenc^j result. 



Theorem 3 There is a sequence of feasible solutions to problem such that its 
w-component converges from above to Ix in L 1 norm and almost uniformly. 

Proof: Follows by the same arguments as Theorem 3 in [6] . □ 



6 LMI relaxations 

In this section we derive finite dimensional semidefinite programming (SDP) or linear ma- 
trix inequality (LMI) relaxations to the infinite dimensional LPs Q and Q and establish 
several convergence results relating these relaxations to the infinite dimensional LPs and 
to the ROA. 

In what follows, Kfc[-] denotes the vector space of real multivariate polynomials of total 
degree less than or equal to k. 

Derivation of the finite dimensional relaxations is standard and the reader is referred to [6J 
Section 5] or to the comprehensive reference [9]; therefore we only highlight the main 
ideas. First of all, since the supports of all measures feasible in Q are compact, these 
measures are determined by their moments, i.e., by integrals of all monomials (which is a 
sequence of real numbers when indexed in, e.g., the canonical monomial basis). Therefore, 
it suffices to restrict the test functions w(x) and v(t,x) in Q to all monomials, reducing 
the linear equality constraints of to linear equality constraints on the moments. Next, 
by the celebrated Putinar Positivstellensatz (see 0QJ]), the constraint that the support 
of a measure is included in a given compact basic semialgebraic set is equivalent to the 
feasibility of an infinite sequence of LMIs involving the so-called moment and localizing 
matrices, which are linear in the coefficients of the moment sequence. By truncating 
the moment sequence and taking only the moments corresponding to monomials of total 
degree less than or equal to 2 k we obtain a necessary condition for this truncated moment 
sequence to be the first part of a moment sequence of a measure with the desired support. 

This procedure leads to the primal SDP relaxation of order k 

p* k = max (y ) 
s.t. 



M k - dT (-g x ,y^)hO 
M k - dT (-g T ,y*)hO 



(9) 

where the notation b stands for positive semidefinite and the minimum is over moment 
sequences (y, yo, y\, y\^ y ) truncated to degree 2k corresponding to measures /i, /i , /i^, 
and jj, . The linear equality constraint captures the two linear equality constraints 



A k (y,yo,yT,yTiyo) = h 








M k (y) b o, 


M k . 


-d Xi (9x 


,2/) b0 


M k (y ) b 0, 
M k (y l T ) b 0, 
M k (y 2 T ) b 0, 


M k . 


-d x (gx: 


,2/o) bo 


M k . 


-d T (gx, 


2/i) b0 
y 2 T ) b0 


M k . 


-d T {gx, 


M k {y ) h 0, 


M k . 


-d x (gx- 


M bo 


M k ^(t(T-t) : y) b0, 


M k . 


-MT- 


-t),y l T ) b0 



10 Please refer to [6] or, e.g., [3] for definitions of the various types of convergence relevant in this 
context. 
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of Q with v(t, x) G M 2 fc[t, x] and G M 2 fc[x] being monomials of total degree less than 
or equal to 2k. The matrices M&(-) are the moment and localizing matrices, following the 
notations of [9j or [6j. In problem ([9]), a linear objective is minimized subject to linear 
equality constraints and LMI constraints; therefore problem pi) is an SDP problem. 



The SDP problem dual to problem ^ turns out to be the sum-of-squares problem 

d* k — inf w ■ I 

s.t. —Cv(t, x) — p(t, x) + qi(t, x)t(T — t) + q 2 {t, x)gx{x) 
w(x) - v(0, x)-l= p {x) + q 0l (x)g x {x) 
v(t, x) = Pti{x) + qri{t, x)t(T -t) + r(x)g x {x) 
v(T, x) = p T2 {x) + Qt 2 ( x )9x(x) - qT^{x)g T {x) 
w(x) = s (x) + si(x)g x (x), 



(10) 



where / is the vector of Lebesgue moments over X indexed in the same basis in which 
the polynomial w(x) with coefficients w is expressed. The minimum is over polynomials 
v(t, x) G M 2 fc[£, x], w(x) G R 2 fc[:c], over the polynomial r(x) and polynomial sum-of-squares 
p(t,x), qi{t,x), q 2 {t,x), q 0l (x), p Tl (x), Pt 2 { x ), Qti(x), Qt 2 {x), q T3 (x), s o{x), s x (x) of 
appropriate degrees. The constraints that polynomials are sum-of-squares can be written 
explicitly as LMI constraints (see, e.g., [9j), and the objective is linear in the coefficients 



of the polynomial w(x); therefore problem (10) can be formulated as an SDP problem 



Theorem 4 There is no duality gap between primal LMI problem f9l) and dual LMI 



problem (10), i.e. p* k = d* k . 



Proof: Follows by the same arguments based on standard SDP duality theory as Theo- 
rem 4 in [6]. □ 

Now we prove our main convergence results. 

Theorem 5 Let Wk G M 2 fc[^] denote the w -component of a solution to the dual LMI 



problem (10) and let Wk{x) = min^ i(^(rr). Then 1 — Wk converges from below to Ix in 



L 1 norm and 1 — Wk converges from below to Ix in L 1 norm and almost uniformly. 

Proof: It follows from Theorem [3] and from the density of polynomials in the space of 
continuous functions on compact sets (for a detailed argument in a similar setting see [6j 
Theorem 5]) that Wk and w~k converge from above to ixg i n Li an d almost uniformly on 
X, respectively. Therefore 1 — Wk and 1 — Wk converge from below to Ix = 1 — Ix c on X 
in the same manner. □ 

The next Corollary follows immediately from Theorem [5j 



Corollary 1 The sequence of infima of LMI problems (10) converges monotonically from 
above to the supremum of the LP problem i.e., d* < d* k+1 < d* k and lim^oo d k = d* . 
Similarly, the sequence of maxima of LMI problems |^) converges monotonically from 



above to the maximum of the LP problem i.e., p* < p%+\ — Pk an< ^ nm fc-^ooPj = P* ■ 

Proof: Follows the proof of Corollary 1 in [6] . Monotone convergence of the dual optima 
d* k follows immediately from Theorem [5] and from the fact that the higher the relaxation 
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order fc, the looser the constraint set of the minimization problem (10). To prove conver- 
gence of the primal maxima observe that from weak SDP duality we have d\ > pi and 
from Theorems [5] and [2] it follows that d* k — >■ d* = p* . In addition, clearly p\ > p* and 
Pk+i — Pk smce the higher the relaxation order k 1 the tighter the constraint set of the 
maximization problem ([9]). Therefore p* k — > p* monotonically from above. □ 

Our last results establishes set-wise convergence of inner approximations to the ROA. 
Theorem 6 Let Wk G M 2 fc[^] denote the w -component of a solution to the dual LMI 



problem (10) and let X 0k :— {x G X : Wk{x) < 1}. Then X 0k C X , 
lim A(X \ X ok ) = and X(X \ U^ =1 X 0fc ) = 0. 



k— >oo 



Proof: Follows the proof of Theorem 6 in [6]. From Lemma [I] we have w^ix) < 1 x G 
Xq for all x G X, and therefore Ix ok < Ix - Since Wk > on X we also have 1 — < Ix ok 
on X and therefore 1 — Wk < Ix ok < Ix on X. From Theorem |5j we have 1 — Wk —> Ix 
in L 1 norm on X. Consequently, 

X(X ) = / I Xo d\ = lim / 1 - w fc c?A < lim / 7x ofc = lim X(X ok ) 

J X k— >oo J j£ k— >oo J j£ k— >oo 

< lim \(utiXoi) = A(U£° =1 X ofc ). 

k— >oo 



But since X ok C X Q for all k, we must have 
lim X(X ok ) = X(X ) ar 

which proves the theorem. □ 



lim A(X ofc ) = X(X ) and A(U^ =1 X *) = A(X ), 

k— >oo 



7 Numerical examples 

In this section we present two numerical examples. The primal problems on measures were 
modeled using Gloptipoly 3 [7] interfaced with the SDP solver SeDuMi [12]; this solver 



also returns the solution to the dual SDP relaxation. In Section |7.3| we then investigate 
how tight low order approximations can be obtained. 

7.1 Univariate cubic dynamics 

Consider the system given by 

x = x(x — 0.5)(x + 0.5), 

the constraint set X = [—1, 1], the final time T — 10 and the target set X T — [—0.3, 0.3]. 
The ROA in this setup is Xq = [—0.5, 0.5]. Polynomial approximations to the complement 
ROA for degrees d G {16,20,24,28} are shown in Figure [l| As expected, the functional 
convergence of the polynomial to the discontinuous indicator function is rather slow. A 
slightly better convergence is observed in the volume error of the sublevel set approxi- 
mation to the ROA documented in Table [T] The relatively slow convergence could be 



significantly improved if a tighter constraint set X was employed; see Section |7.3| be- 
low. Alternative polynomial bases (e.g. Chebyshev polynomials) would also allow tighter 
higher order approximations; see [8j for more details. 
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Figure 1: Univariate cubic dynamics polynomial approximations (solid line) to the 
complement ROA indicator function Ix c = /[_i,-o.5] + ^[0.5,1] (dashed line) for degrees 
d G {16,20,24,28}. 



7.2 Van der Pol oscillator 



As a second example consider a scalec 11 version of the uncontrolled reversed-time Van der Pol 
oscillator given by 

X\ — — 2.0x 2 , 

x 2 = 0.80xi + 10(^ - 0.21)^2. 

We take T = 1 and X T = {x G R n : ||x|| 2 < 0.50} and X := {x G R n : \\x\\ < 1.1}. The 
ROA is bounded, having the characteristic Van der Pol shape. Plots of polynomial sublevel 
set approximations of degrees d G {9, 12, 15, 18} are shown in Figure |3[ We observe a 
relatively fast convergence to the ROA, which is also documented by the relative volume 
errors reported in Table [2j Figure [2] then shows a degree 18 polynomial approximation to 
the indicator function of the complement ROA. 



11 The coefficients were chosen so that the ROA fits within the box [—1,1]' 
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Table 1: Univariate cubic dynamics - relative error of the inner approximations to the ROA 
Xq = [—0.5, 0.5] as a function of the approximating polynomial degree. 



degree 


16 


20 


24 


28 


error 


11.4% 


6.4% 


4.84% 


4.54% 




Figure 2: Van der Pol oscillator degree 18 polynomial approximation to the indicator 
function of the complement ROA. 



Table 2: Van der Pol oscillator - relative error of the inner approximation to the ROA Xq as a 
function of the approximating polynomial degree. 



degree 


9 


12 15 


18 


error 


18.3% 


8.4% 3.8% 


3.1% 



7.3 Low order approximations 

In the examples above, relatively high order polynomials had to be used to obtain tight 
approximations, which can limit subsequent applicability of the approximations. There 
are several ways to obtain low order approximations of similar quality. First of all, since 
the integral of a polynomial w is minimized over the constraint set X 1 it is desirable that X 
be a good outer approximation of the ROA. Of course, selecting X is possible only if it is 
an artificially specified outer approximation of the ROA, not a constraint set coming from 



physical requirements on the system. More importantly, notice that in problem (10) the 
system dynamics enters the constraints on the polynomial v(t, x), whereas the polynomial 
w(x) is only upper-bounding v(t,x) + 1 for t — 0. Since the inner approximations are 
given by sublevel sets of w, it is possible and plausible to choose different degrees of w and 
v - low for w and higher for v. Both techniques are illustrated in Figure [4j in Figure [4] (a) 
we consider the univariate cubic dynamics and we both shrink the constraint set X and 
choose low order w while keeping v of higher order. In Figure [4] (b) we consider the Van 
Der Pol oscillator, keeping the constraint set X unchanged and only selecting low order w. 
The inner approximations obtained are indeed significantly tighter for the given degrees 
(compare with Figures [I] and |3|) . 
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-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



X 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



X X 

Figure 3: Van der Pol oscillator polynomial inner approximations (light gray) to the 
ROA (dark gray) for degrees d G {9, 12, 15, 18}. 




~_1 _0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



(a) Univariate cubic dynamics - constraint set (b) Van der Pol oscillator - deg w — 8 (deg v — 
X — [ — 0.7,0.7], degu> = 6 (degv = 16). Volume 18). Volume approximation error 5.46%. 
approximation error 2.25%. 

Figure 4: Low order approximations to the ROA. Left: tighter constraint set X and low 
order w (compare with Figure flT) . Right: low order w only (compare with Figure p|). 
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8 Conclusion 



This paper presented an infinite dimensional convex characterization of the region of at- 
traction (ROA) for uncontrolled polynomial systems, following the approach initiated in 
our previous work [6j . Finite dimensional dual relaxations yield a converging sequence of 
inner approximations to the ROA, thereby complementing the outer approximations of 
[6] . One of the virtues of the approach is its conceptual simplicity the resulting approx- 
imation is the outcome of a single SDP or LMI problem with no free parameters except 
for the relaxation order. The approximations itself are also simple, given by sublevel sets 
of polynomials of predefined degrees. 

Nevertheless, this approach does not escape the curse of dimensionality indeed, whereas 
the number of variables of the LMI relaxations grows polynomially with the relaxation 
order, this number grows exponentially with the state dimension. Tailored structure- 
exploiting SDP solvers could enable this approach to reach higher dimensions. In addition, 
a different choice of basis functions (e.g., Chebyshev polynomials rather than monomials) 
would improve numerical conditioning of the LMIs, allowing higher oder relaxations to 
be computed. 

Future research directions include inner approximations in a controlled setting and the 
related problem of robust region of attraction / reachable set computation with either 
unknown but constant uncertainty or a time- varying disturbance. The cases of asymptotic 
region of attraction and maximum (controlled) positively invariant set computation are 
amenable to similar tools. 
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